/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#if !defined(_2D_Representation_)
#define _2D_Representation_

#include <src/data/ParticleSet.h>
#include <src/data/OS_Multi_Array.h>
#include <src/math/Operators/particle_density.h>

template<class data, std::size_t = data::Dim, typename Family_Operator = Null>
struct representation_2 {
	typedef particle<typename data::base_value_type>				p_array_1;
	typedef particles_list<p_array_1>								p_array_2;
	p_array_2 operator()(const data & from) {OS_DYNAMIC_CHECK(false,"representation_2 : Overload this function for your own set of data");}; // to implement
};

// we assume at present time that the data are computed in [0,1]
template<class data>
struct representation_2<data,1,Finite_Diff_Operator> {
	typedef particle<typename data::base_value_type>				p_array_1;
	typedef particles_list<p_array_1>								p_array_2;
	p_array_2 operator()(const data & from) {
		OS_STATIC_CHECK(data::Dim == 1); // "representation_2 is not implemented for this type of data";
		p_array_2 temp(from.size(),2);
		for (OS_size i = 0; i < temp.size(); ++i) {
			temp[i][1] = from[i];
			temp[i][0] = (i+0.5) / temp.size();
		}
		return temp;
	};
};

template<class data>
struct representation_2<data,1,Particle_Operator> {
		typedef particle<typename data::base_value_type>				p_array_1;
		typedef particles_list<p_array_1>								p_array_2;
		p_array_2 operator()(const data & from) {

			data density_ = particle_density<data>()(from);

			p_array_2 temp(from.size(),2);
			for (OS_size i = 0; i < temp.size(); ++i) {
				temp[i][0] = from[i];
				temp[i][1] = density_[i];
			};
			return temp;
	};
};

template<class data>
struct representation_2<data,2,Particle_Operator> {
		typedef particle<typename data::base_value_type>				p_array_1;
		typedef particles_list<p_array_1>								p_array_2;
		p_array_2 operator()(const data & from) {
			p_array_2 temp(from.size(2),2);
			for (OS_size i = 0; i < temp.size(); ++i) {
				temp[i][1] = from[i][1];
				temp[i][0] = from[i][0];
			};
			return temp;
	};
};
template<class data>
struct representation_2<data,2,Finite_Diff_Operator> {
		typedef particle<typename data::base_value_type>				p_array_1;
		typedef particles_list<p_array_1>								p_array_2;
		p_array_2 operator()(const data & from) {
			p_array_2 temp(from.size(1),2);
			for (OS_size i = 0; i < from.size(1); ++i) {
				temp[i][0] = (i+0.5) / from.size(1);
/*				data::base_value_type first		= from[0][i];
				data::base_value_type second	= from[1][i];
				temp[i][1] = std::sqrt(first*first + second*second);
*/
				temp[i][1] = from[0][i];
			};
			return temp;
	};
};


template<class data, std::size_t = data::Dim, typename Family_Operator = Null>
struct fd_to_p__representation_2 {
	typedef particle<typename data::base_value_type>				p_array_1;
	typedef particles_list<p_array_1>								p_array_2;
	p_array_2 operator()(const data & from) {OS_STATIC_CHECK(false);}; // to implement
};

template<class data>
struct fd_to_p__representation_2<data,1,Finite_Diff_Operator> {
			typedef particle<typename data::base_value_type>				p_array_1;
			typedef particles_list<p_array_1>								p_array_2;
			p_array_2 operator()(const data & from) {

			data average_ = particle_average<data>()(from);
			data density_ = ::inv_particle_density<data>()(from);

			p_array_2 temp(average_.size(),2);
			for (OS_size i = 0; i < temp.size(); ++i) {
				temp[i][0] = density_[i];
				temp[i][1] = average_[i];
			};
			return temp;
	};
};


template<class data, std::size_t = data::Dim, typename Family_Operator = Null>
struct representation_3_common {
	typedef particle<typename data::base_value_type>				p_array_1;
	typedef particles_list<p_array_1>								p_array_2;
	virtual p_array_2 operator()(const data & from) {
		OS_DYNAMIC_CHECK(false,"representation_3_common::operator() implment this function");
		return p_array_2();
	}; // to implement
};

template<class data, std::size_t = data::Dim, typename Family_Operator = Null>
struct representation_3 {
	OS_STATIC_CHECK(false); // to implement
};

template<class data>
struct representation_3<data,2,Finite_Diff_Operator> : public representation_3_common<data,2,Finite_Diff_Operator>{
		virtual p_array_2 operator()(const data & from) {
			OS_size size_1 = from.size(1);
			OS_size size_2 = from.size(2);
			OS_double alpha = .1;
			OS_double coef1 = std::sin(alpha);
			OS_double coef2 = std::cos(alpha);
			p_array_2 temp(size_1*size_2,2);
			for (OS_size i = 0; i < size_1; ++i) {
				for (OS_size j = 0; j < size_2; ++j) {
					double breakpoint = from[i][j];
					temp[i+j*size_1][0] = ( (i+.5) / size_1 - j*coef1/ size_1);
					temp[i+j*size_1][1] = (j+.5) * coef1/ size_2 + breakpoint;
				};
			};
			return temp;
	};
};





#endif
